library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(readxl)
library(leaflet)
options(
tigris_class = "sf"
)
Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")
Mapping population of ages 65 and older for the Bay Area
# If not already done, obtain the census data
# acs_vars <-
# listCensusMetadata(
# name = "2018/acs/acs5",
# type = "variables"
# )
# # Save an .rds version on my computer
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
# saveRDS(acs_vars, file = "censusData2018_acs_acs5.rds")
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# obtain the saved census data
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# define Bay Area county names
bay_county_names <-
c(
"Alameda",
"Contra Costa",
"Marin",
"Napa",
"San Francisco",
"San Mateo",
"Santa Clara",
"Solano",
"Sonoma"
)
# get Bay Area tracts
bay_tracts <-
bay_county_names %>%
map(function(x){
tracts("CA",x,progress_bar=F) %>%
pull(GEOID)
}) %>% unlist()
# get sex by age for Bay Area counties, modified from the rBasics.Rmd code
bayArea_sexbyage <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "tract:*",
regionin = "state:06",
vars = "group(B01001)"
) %>%
mutate(
tractval =
paste0(state,county,tract)
) %>%
filter(tractval %in% bay_tracts) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(
key = "variable",
value = "estimate",
-tractval
) %>%
mutate(
label = acs_vars$label[match(variable,acs_vars$name)]
) %>%
select(-variable) %>%
separate(
label,
into = c(NA,NA,"sex","age"),
sep = "!!"
)
# Now find the population that is age 65 and up - also modified from rbasics.Rmd
bayArea_elderly <-
bayArea_sexbyage %>%
filter(!is.na(age)) %>%
mutate(
elderly =
ifelse(
age %in% c(
"65 and 66 years",
"67 to 69 years",
"70 to 74 years",
"75 to 79 years",
"80 to 84 years",
"85 years and over"
),
estimate,
NA
)
) %>%
group_by(tractval) %>%
summarize(
elderly = sum(elderly, na.rm = T),
total_population = sum(estimate, na.rm = T)
) %>%
mutate(
percent_elderly = elderly/total_population
)
# get Bay Area tracts objects
bay_tracts_obj <- NULL
for (i in 1:length(bay_county_names)) {
bay_tracts_obj <- bay_tracts_obj %>% rbind(tracts("CA",bay_county_names[i],progress_bar=F))
}
# plot percent elderly
bayArea_elderly %>%
left_join(bay_tracts_obj %>% dplyr::select(GEOID), by = c("tractval" = "GEOID")) %>%
st_as_sf() %>%
filter(!is.na(percent_elderly)) %>%
mapview(zcol = "percent_elderly")